Phase diagram of the rotating two-component Fermi gas including vortices 
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We determine the conditions under which superfluidity with and without quantized vortices appears in a 
weakly interacting two-component atomic Fermi gas that is trapped in a rotating cylindrical symmetric harmonic 
potential. We compute the phase diagram as a function of rotation frequency, scattering length, temperature, 
total number of trapped atoms, and population imbalance. 



Introduction. Superfiuids are fluids that can flow with 
hardly any friction. If a superfluid is put into a rotating con- 
tainer, vortices carrying quantized circulation can be formed. 
These vortices are the hallmark of superfluidity and have 
been observed experimentally in Helium-4, in atomic Bose- 
Einstein condensates and in a two-component atomic Fermi 
gas made out of 6 Li atoms E[2). 

In this Letter we will focus on vortex formation in an equi- 
librated weakly interacting two-component Fermi gas that is 
trapped in a rotating cylindrical symmetric harmonic poten- 
tial. The two components of this gas consist of atoms in differ- 
ent hyperfine states, and will be labeled by f , 4- We will con- 
sider a situation in which both components have equal mass 
M. In experiment, one can vary the strength of the interaction 
between the components using the Feshbach resonance. Fur- 
thermore one can control the rotation frequency £2, tempera- 
ture T, total number of atoms, and the relative difference be- 
tween the number of atoms in each component P (population 
imbalance) |Q][2j. When the interaction is tuned to be attrac- 
tive, the gas is cooled to low enough T, and Q. and P are made 
small enough, the components will form Cooper pairs result- 
ing in a Bardeen-Cooper-Schrieffer (BCS) superfluid 0. 

The aim of this Letter is to determine in which region of 
the parameter space that is accessible to experiment vortices 
will be formed. So far only very limited information about 
this region is available. Only for weak interactions, T = 
and P = 0, the lower critical Q. has been obtained theoreti- 
cally. This frequency was estimated in [4| using a Ginzburg- 
Landau approach and computed by the author through solving 
the Bogoliubov-de Gennes (BdG) equation [5]. Experimen- 
tally only the critical P for vortex formation at one specific 
Q. has been determined for different scattering lengths in the 
strongly interacting regime (2)- 

We will consider the following trapping potential: U(x) - 
jMaj 2 p 2 , where we have introduced cylindrical coordinates, 
i.e. x - (p cos <p,p sin <f>,z). We will allow the trap to rotate 
with constant angular frequency Q. in the x-y plane. This trap- 
ping potential implies harmonic confinement with frequency 
a) in the x-y plane and infinite extension in the z -direction. In 
an experiment this situation can be approached by choosing 
the trapping frequency in the z-direction (a> z ) much smaller 
than u>. The characteristic length scale of the potential is 
the harmonic oscillator length A = y/h/Maj. In the setup 
of Ref. |2i| the radial trapping frequency was taken to be 



lu/(2jt) =110 Hz. In that case A ~ 3.9 pm and the char- 
acteristic energy scale hu>/k B ~ 5.3 nK. 

The order parameter for superfluidity is the pairing field 
A(x). For a single vortex that is located at the center of the 
trap it has the following form: A(x) = A{p)exp(ik<p), with 
k the winding number of the vortex. The k = case corre- 
sponds to a superfluid without vortices. We will assume that 
A(p) e R. There is no superfluidity at the core of a vortex, 
hence A(p = 0) = for k + 0. If Q. = 0, T = 0, and P = 0, the 
whole system of atoms forms a vortex-free superfluid. In that 
situation the vortex is metastable. This is because condensa- 
tion energy is lost near the vortex core, the atoms traveling 
around the vortex core have an average velocity resulting in a 
kinetic energy cost, and the system has to expand to compen- 
sate for the density depletion at the vortex core. This expan- 
sion costs energy because of the attractive interaction. 

Starting from a situation in which Q = 0, let us now imag- 
ine increasing Q. If both T = and P = 0, the whole system 
stays superfluid up to a critical frequency Q/,. Above Q;, part 
of the system will turn into a normal gas due to breaking of 
Cooper pairs (6). If either T + or P + there are un- 
paired atoms present for any nonzero D.. The unpaired atoms 
are predominantly located in the outer regions of the system. 
They rotate like a rigid body and acquire therefore rotational 
energy. A vortex is also a source of rotational energy since 
it induces angular momentum in the system. If the rotational 
energy gains overcome the costs, a single vortex will be pre- 
ferred above a frequency Q;. Because of symmetry and energy 
arguments, this vortex has unit winding number (k = 1) and 
is located at the center of the trap (p = 0). By further increas- 
ing Q more vortices can be created resulting in a vortex lattice 
EH). At the same time increasing Q. will shrink the size of 
the superfluid region [9], while the system as a whole will ex- 
pand due to the centrifugal force. At some point the superfluid 
region will be so small that it can only support a single vortex 
at p = 0. By further increasing O this vortex will disappear at 
a upper critical rotation frequency Q„. Then at an even larger 
rotation frequency Q, superfluidity will vanish completely via 
a second order transition [7]. The upper critical frequency for 
superfluidity Q, has been computed in Refs. |9, 10 1 for a bal- 
anced gas using the local density approximation (see also ifTTIl 
for related analyses of the phase boundary of superfluidity). 
In this Letter we will obtain Cl s by solving the BdG equation. 
Finally, above Q. — a> the system will be torn apart. 
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As pointed out above, the first vortex that appears and the 
last vortex that disappears when increasing Q has k — 1 and 
is located at p = 0. Moreover, in experiment it has been ob- 
served that the number of vortices goes continuously to zero 
when increasing P at fixed Q [2|. Therefore, we can map 
out the entire region in which one or more vortices are pre- 
ferred thermodynamically by determining the conditions un- 
der which the k = 1 situation has lower Helmholtz free energy 
than the k — situation. For studies of real-time dynamics of 
vortex formation we refer to Refs. Il8l [T21 . To obtain Q s we 
will determine the point at which A(p = 0) vanishes for k — 0. 

Setup. We will now briefly discuss the details of the calcu- 
lation, a more extensive discussion can be found in our earlier 
work 0. To obtain the Helmholtz free energy one first needs 
to compute the pairing field A(jc) = A(p) exp(ik<p) and the den- 
sity profiles pn(x) = Pfiip) for each component separately. 
For a weakly interacting gas these can be found by solving 
the BdG equation self-consistently: 
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Here 7^ contains the single-particle Hamiltonian, the 
Hartree self-energy, and the chemical potential p.^. Explic- 
itly, it reads Wu(fl) = ^ + iJWw 2 p 2 - OL z - Pj^+S^uip), 
where the z-component of the angular momentum is given by 
L z = -iHd/d(f), and g = 4nah 2 /M is the coupling constant with 
a the s-wave scattering length. The wavefunctions u,(x) and 
v,(jc) have to be normalized as Jd 3 x [\m(x)\ 2 + |v;(*)| 2 ] = 1. 
The number densities are given by tif(x) = 2; f(Ei)\ui(x)\ 2 
and ni (x) = Zif(-EM(x)\ 2 , where f(E) = [exp(££) + l]" 1 
with j6 = l/(ksT). The pairing field follows from the regular 
(reg) part of the anomalous propagator in the following way 
A(x) = gG T ^(x,x) where G n (*,x') = I li f(E i )u i (x)v*(x'). To 
obtain this regular part we have used a method |5| based on 
the procedures discussed in Refs. |fT3l[T4ll . 

We will consider a fixed number of atoms per unit length 
A in the z-direction, and denote this number by A/|,j, = 
A/Lj d?xn^i(x). Here L is the length of the system in 
the z-direction which we will take to be infinite. We will 
write N - N-\ + N[ for the total number of atoms per unit 
length. The population imbalance or polarization is defined as 
P = (Nf—N{)/N. The chemical potentials p.^ will be solved 
for such that the required N^.i is obtained. 

To solve the BdG equation numerically, we have discretized 
the radial part of the wavefunctions on a Lagrange mesh 
|[T5l based on Maxwell polynomials [5 |. Typically we could 
reach a relative accuracy of order 1CT 3 with about 64 to 96 
mesh points for N = 1000. The angular and z-dependence 
of the wavefunctions were treated exactly. Integration over 
the z-momentum was performed using the adaptive Simp- 
son method. To solve for self-consistency we have used the 
Newton-Broyden rootfinding method. The full details of our 
numerical procedure are explained in J3]- Examples of pair- 
ing field and density profiles with and without a vortex can be 
found in Refs. [T3] [16] [TTJ . In particular we have obtained 




FIG. 1: Phase diagram: a-Q. plane, for T ~ 0, Af = 1000 and P = 0. 
The dotted line indicates Q b . Here 0.8 < l/(k F0 \a\) < 2.6. 



excellent agreement with Ref. ifTTll . in which a vortex profile 
for an imbalanced gas that is also trapped in a cylindrical sym- 
metric harmonic potential is presented. 

Once the pairing field and the density profiles have been 
obtained, the Helmholtz free energy per unit length (which is 
ultraviolet finite) can be computed in the following way 



i L " 

- ^ J d 3 x [G n (x, xYA(x) + gn t (x) ni (*)] + | J] e u (2) 



where e,- are the eigenvalues of the Hartree-Fock Hamilto- 
nian H H f = VH^D. = 0) + tii(Cl = 0)] 12. We have computed 
AT - Tk=\ —fk=o for several values of the external parameters 
and obtained the phase boundary by determining the point at 
which AT = 0. Typically we could locate this boundary with 
a relative accuracy of about 10~ 2 to 10~ 3 . 

Results. We will now present several phase diagrams from 
which it can be seen for which values of the rotation frequency 
Q., scattering length a, temperature T, total number of atoms 
per unit length N, and imbalance P, superfluidity with one or 
more vortices will be formed (light gray area). Also we will 
indicate in these phase diagrams when the system exhibits su- 
perfluidity without any vortex (dark gray area). As a measure 
of the interaction strength we will use l/\kfoa\, where kfo de- 
notes the Fermi momentum of the superfluid at T — 0, P = 0, 
£2 = and p = 0. 

In Ref. J5], we have computed the critical frequency for 
unpairing (Q/,) and the lower critical frequency for vortex for- 
mation (Q/) as a function of scattering length for T — 0, P — 
and N = 1000. In Fig. [T] we display the full phase diagram 
at T ~ 0. Here we write T ~ to indicate that the vortex 
phase boundaries and the unpairing transition were computed 
at T — exactly, whereas Q s was computed at T — Q.QITioj. 
We have used this very small T to ensure that A(p = 0) ap- 
proaches zero continuously when increasing Q, so that we 
could determine Q s . Exactly at T = A(p = 0) does not 
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FIG. 2: Phase diagrams: N-Q. plane, for T ~ 0, P = 0, and 
1/feoM) = 1.2 (left) and 1/feoM) = 0.8 (right). 
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FIG. 3: Phase diagrams: Q-T plane, for N = 1000, P = and a = 
-0.10A (left, 1/feoM) = 1.2) and a = -0.13/1 (right, l/(k F0 \a\) = 
0.8). The dot indicates the unpairing quantum critical point. 



seem to vanish when increasing £2, although it does become 
very small. 

The successive transitions that one encounters when in- 
creasing £2 were already described in the introduction and can 
be clearly seen in the diagram. The first transition is the un- 
pairing transition occurring at £2 = £1/,. It is of second order 
and turns into a cross-over for T > 0. Therefore it is a quan- 
tum phase transition and at £2/, the system resides at a quantum 
critical point. If \a\ increases, it will become more difficult 
to break the Cooper pairs, hence Q.t, grows in that case. For 
\a\ > 0.085/1 vortices will be formed for £2/ < £2 < £2„. The 
structure of £2/ is a result of the interplay of two effects J3]. 
The first is that the energy cost of a vortex at Q = increases 
when increasing \a\. That naturally leads to a larger £2/. The 
second effect is that unpairing becomes easier for smaller \a\, 
hence more rotational energy can be acquired by the k = su- 
perfluid. This leads to an increase of £2/ for weak interactions 
and is the reason that for small \a\ no vortices will be formed 
for any £2. The size of the superfluid region shrinks above £2/, 
when increasing Q. If \a\ grows at fixed Q, it becomes more 
difficult to destroy superfluidity by rotation. This results in a 
larger superfluid region, so that a vortex can fit more easily. 
For these reasons both £2/ and the upper critical frequency for 
superfluidity £2, grow with increasing \a\. The behavior of £2 S 
is qualitatively in agreement with the results of Refs. ll9l [T0l . 

This phase diagram will be modified quantitatively when 
changing N. However, the qualitative structure will remain 
the same if N is large enough. To illustrate this, we display 
in Fig. [2] the phase diagram in the W-Q plane for two differ- 
ent interaction strengths. A larger N implies a larger system, 
and hence for a given £2 the atoms at the boundaries have a 
larger velocity. This makes pairing more difficult, leading to a 
smaller Q.t, as can be seen in Fig. [2] To explain the behavior of 
£2/, we can use that the rotational energy gain of the vortex is 
proportional to the angular momentum which is proportional 
to N. The energy costs of the vortex grow much slower when 
increasing N. Hence a larger N leads to a smaller £2/. Below 
a certain N no vortices will be formed for any Q. The upper 
critical frequencies £2„ and £2, increase if N grows. 

If T is increased the window in Q in which vortices are 
formed narrows. This can be seen from Fig. [3] in which we 
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FIG. 4: As in Fig.[3] but now in Q-P plane for T = 0. 

display the phase diagram in the Q.-T plane for N = 1000, 
P — and two different scattering lengths. Above a certain 
T, vortices will not be formed for any Q. while the system is 
still partly superfluid. It is the easiest to make vortices at in- 
termediate £2, since in that case the window in T is the largest. 
This Q.-T phase diagram is qualitatively very similar to that of 
a rotating Bose-Einstein condensate [ 18 1. 

In Fig. [4] we display the phase diagram in the D.-P plane for 
T = and two different scattering lengths. The larger the P, 
the narrower the window in £2 in which vortices are formed 
becomes. Above a certain P part of the system can still be 
superfluid, but vortices will not be formed for any Q. 

From Figs. [3] and [4] it can be seen that the critical T and 
critical P for vortex formation and superfluidity grow when 
increasing the interaction strength. When Q is increased, Q s 
always decreases. Especially for large interaction strengths 
Q„ and £2 S lie very close each other. 

In the experiment described in Ref. a strongly inter- 
acting two-component Fermi gas made out of ~ 7 x 10 6 
atoms was trapped in a cigar-shaped potential with a fre- 
quency ratio of u> z /(l> = 23/110. Using the Thomas-Fermi 
approximation at T — we estimate that in this experiment 
N(z = 0) ~ 1 x 10 5 . The atoms were stirred with a laser 
at a frequency of Q = (70/1 10)w ~ 0.636w. After stirring, 
one waited until the system had equilibrated and measured the 
number of vortices as a function of P. In this way an upper 
critical imbalance for vortex formation (P c ) could be deter- 
mined. Although this situation is not completely equivalent to 
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FIG. 5: Phase diagram: P-T plane, for N = 1000, Q. = 0.636w, and 
a = -0.13/1. Here l/feuM) = 0.8. 



have obtained are quantitatively reliable and are in principle 
directly comparable to a possible future experimental deter- 
mination. Our analysis can be extended to more complicated 
systems, like Fermi gases with /j-wave pairing, Fermi gases 
with more than two components, and Fermi gases in which 
the two components have unequal mass. This will be useful 
for the experimental search for superfluidity in such systems. 
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FIG. 6: Phase diagrams: N-P plane, for T = 0, Si = 0.636ai, and 
1/feoN) = 1-0 (left) and l/(k P0 \a\) = 0.8 (right). 

our setup it is nevertheless interesting to make a comparison 
to this experiment. Therefore we display in Fig. [5] the phase 
diagram in the P-T plane for Q = 0.636o>, N — 1000, and 
l/(kfo\a\) - 0.8. In qualitative agreement with the experiment 
it can be seen that P c is weakly dependent on T for small tem- 
peratures. We find that for P > P c there is a small window 
in which part of the system is in the superfluid phase without 
any vortex, as in the experiment. 

In Fig.|6]we display the phase diagram in the N-P plane for 
T — 0, Q = 0.636w and two different interaction strengths. As 
expected and seen in experiment, weaker interactions imply 
a lower P c . In the experiment it was found that P c ~ 0.1 
at l/(kfo\a\) = 0.5 0. The values of P c we have obtained 
seem to be large compared to this value, because we consider 
weaker interactions and our N(z = 0) is much smaller. This 
quantitative difference could have been caused by the different 
stirring method or by the shape of the trapping potential. Also 
it could signal that one should take into account beyond the 
mean field corrections at l/(k F o\a\) = 0.8. 

Conclusions. For the first time we have unveiled the full 
phase diagram of a trapped, rotating, and weakly-interacting 
two-component Fermi gas including vortices. We have made 
detailed predictions for the conditions under which superflu- 
idity with and without vortices is formed as a function of ro- 
tation frequency, scattering length, temperature, number of 
atoms and population imbalance. The phase diagrams we 
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